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We  present  and  use  an  algorithm  based  on  convex  conic  optimization  to  design  two-dimensional  photonic 
crystals  with  large  absolute  band  gaps.  Among  several  illustrations  we  show  that  it  is  possible  to  design  photonic 
crystals  which  exhibit  multiple  absolute  band  gaps  for  the  combined  transverse  electric  and  magnetic  modes.  The 
optimized  crystals  show  complicated  patterns  which  are  far  different  from  existing  photonic  crystal  designs.  We 
employ  subspace  approximation  and  mesh  adaptivity  to  enhance  computational  efficiency.  For  some  examples 
involving  two  band  gaps,  we  demonstrate  the  tradeoff  frontier  between  two  different  absolute  band  gaps. 
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I.  INTRODUCTION 

Photonic  crystals  are  periodic  structures  created  from 
the  arrangement  of  low  and  high  index  materials.  They 
are  designed  to  affect  the  motion  of  light  by  prohibiting 
the  propagation  of  electromagnetic  waves  in  all  directions 
within  certain  frequency  ranges  known  as  absolute  band  gaps 
(ABGs).  A  complete  band  gap  (CBG)  refers  to  the  case  when 
the  the  ABG  is  independent  of  the  polarization  of  the  wave. 
Photonic  band  gap  structures  have  proven  very  important  as 
an  enabling  tool  for  the  design  and  fabrication  of  many  novel 
devices  including  frequency  filters,  waveguides,  switches,  and 
optical  buffers  [1-4].  Therefore  the  ability  to  design  materials 
which  have  prescribed  band  gap  diagrams  is  very  important 
from  a  practical  perspective. 

It  is  well  known  that  a  low-index-hole  (e.g.,  air-hole) 
two-dimensional  (2D)  photonic  crystal  has  larger  ABGs  in 
transverse  electric  (TE)  modes,  while  a  high-index-hole  2D 
photonic  crystal  has  larger  ABGs  in  transverse  magnetic 
(TM)  modes.  A  detailed  explanation  for  this  behavior  can 
be  found  in  [5].  This  observation  has  been  extensively  used  to 
create  a  wide  variety  of  photonic  structures  [6,7].  However, 
the  structures  created  using  parametric  studies  combined 
with  physical  reasoning  are  in  general  not  optimal  and 
larger  band  gaps  can  often  be  obtained  when  using  formal 
topology  optimization  methods.  Previous  topology  approaches 
include  gradient-based  approaches  [8-10]  and  evolutionary 
methods  [11,12].  Although  these  methods  have  been  used  to 
produce  useful  designs,  the  band  gap  optimization  problem 
is  a  difficult  nonconvex  optimization  problem,  and  first-order 
(and  other  gradient-based)  methods  suffer  from  low  regularity 
and  nondifferentiability  due  to  the  presence  of  eigenvalue 
multiplicities. 

Recently,  the  work  of  [13]  has  used  extensive  numerical 
optimization  to  produce  designs  of  photonic  structures  with 
very  large  ABGs  for  the  first  15  bands.  In  particular,  the  authors 
therein  propose  a  simple  geometric  scheme  that  provides 
structures  with  very  large  gaps  between  any  two  bands.  They 
conjecture  that  the  globally  optimal  structure  for  TM  modes  is 
the  triangular  distribution  of  circular  rods  and  that  the  globally 
optimal  structure  for  TE  modes  is  a  triangular  distribution  of 


hexagonal  (instead  of  the  commonly  used  circular)  holes,  i.e., 
the  honeycomb  structure.  Interestingly,  corroborative  evidence 
for  this  line  of  conjecture  can  be  found  in  the  current  paper’s 
antecedent  work  [14],  which  used  semidefinite  programming 
and  subspace  methods  as  the  enabling  tools  for  more  general 
(though  still  nonconvex)  design  optimization.  However,  it 
should  be  noted  that  the  geometric  scheme  in  [13]  is  not 
applicable  to  the  combined  TE  and  TM  (CBG)  gaps,  and  in 
fact,  using  extensive  optimization,  no  systematic  approach  to 
CBG  design  was  found  in  [13]. 

In  this  paper,  we  aim  to  optimally  design  photonic  crystals 
which  possess  ABGs  in  several  frequency  bands  (i.e.,  multiple 
ABGs)  and  in  combined  TE  and  TM  modes  (i.e.,  multiple 
CBGs).  A  photonic  crystal  structure  exhibiting  multiple  CBGs 
is  of  considerable  interest  because  it  enables  novel  photonic 
devices  to  operate  with  a  wider  range  of  forbidden  frequencies 
and  in  both  TE  and  TM  modes.  For  instance,  a  material 
for  which  propagation  is  forbidden  at  integer  multiples  of  a 
fundamental  frequency  could  prove  useful  for  the  design  of 
resonant  cavities.  Our  optimization  algorithm  is  an  extension 
of  the  work  in  [14],  which  uses  convex  conic  (semidefinite)  op¬ 
timization  as  a  subroutine  in  a  broader  nonconvex  optimization 
scheme  that  aims  to  compute  optimal  solutions.  We  consider 
both  square  lattice  and  hexagonal  lattice  arrangements.  Our 
algorithm  computes  crystal  design  patterns  which  may  be  very 
different  from  existing  photonic  crystals,  yet  many  of  which  are 
simple  enough  to  be  fabricated  using  current  state-of-the-art 
technology.  Moreover,  the  relative  efficiency  of  our  algorithm 
allows  extra  exploration  of  the  design  space  and  thus  increases 
the  likelihood  that  larger  local  (and  possibly  global)  optimal 
solutions  are  computed.  For  some  examples  involving  two 
band  gaps,  we  estimate  the  tradeoff  frontier  to  demonstrate  the 
tradeoff  in  design  between  two  different  absolute  band  gaps. 

Extensive  analysis  has  generally  revealed  the  dictum  that 
“the  TM  band  gaps  are  favored  in  lattices  of  isolated  high-6 
regions,  and  TE  band  gaps  are  favored  in  connected  lattices” 
[5],  and  indeed  this  has  been  validated  in  [10,13,14]  for 
single-gap  photonic  crystals.  The  computation  of  multiple 
ABG  designs  herein  further  validates  this  statement  for  a  large 
number  of  examples  except  for  the  TE  case  where  we  obtain 
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nonconnected  structures  supporting  multiple  band  gaps.  This 
observation  raises  a  dilemma  for  the  polarization-independent 
CBGs  as  it  seems  that  both  the  isolation  of  high-6  and  the 
connectedness  of  the  structure  can  be  satisfied  only  in  a 
triangular  lattice  arrangement  [5].  The  recent  work  of  [15] 
has  investigated  both  absolute  and  complete  band  gaps  for 
two-dimensional  semiconductor-dielectric  photonic  crystals 
as  a  function  of  the  filling  fraction.  By  considering  photonic 
crystals  made  up  of  doped-GaAs  cylinders  (6r  =  12.8)  em¬ 
bedded  in  a  vacuum  background  (ea  =  1),  the  authors  therein 
found  that  the  semiconductor-dielectric  structure  supports 
multiple  ABGs  at  the  filling  fraction  of  47.5%  for  a  square 
lattice  and  at  the  filling  fraction  of  52.5%  for  a  hexagonal 
lattice.  However,  they  could  find  only  a  single  CBG  for  a 
specific  filling  fraction  and  the  magnitude  of  the  gap-midgap 
ratio  of  approximately  1%  is  very  small.  In  this  paper,  we 
demonstrate  that  for  the  same  materials,  our  optimization 
algorithm  yields  not  one,  but  two  complete  band  gaps  with 
significantly  larger  gap-midgap  ratios,  on  the  order  of  6%  or 
more. 

The  paper  is  organized  as  follows.  In  Sec.  II,  we  present  the 
photonic  crystal  design  optimization  problem  and  present  our 
algorithm  for  computing  large  band  gap  designs.  In  Sec.  Ill,  we 
present  a  mesh  adaptivity  methodology  as  well  as  a  summary 
of  our  computational  scheme.  In  Sec.  IV,  we  present  selected 
results,  including  adaptive  computational  mesh,  crystal  design 
with  multiple  band  gaps  in  TE,  TM,  and  combined  TEM 
polarizations,  and  in  square  as  well  as  hexagonal  lattices. 

II.  PHOTONIC  CRYSTAL  DESIGN  VIA  CONSTRAINED 
OPTIMIZATION 

We  first  review  basic  terminology,  notation,  and  modeling 
framework  from  [14].  A  two-dimensional  photonic  crystal  is 
characterized  by  having  a  dielectric  function  e(r)  which  is 
periodic  in  the  xy  plane  and  constant  in  the  z  direction,  i.e., 
e(r)  =  6(r  +  Rj),  where  Rj,  ( d  =  1,2),  are  primitive  lattice 
vectors  depending  on  the  periodicity  length  a  of  the  lattice.  The 
electromagnetic  modes  of  a  two-dimensional  crystal  can  be 
classified  into  two  different  polarizations:  transverse  magnetic 
(TM)  (electric  field  in  the  z  direction)  and  transverse  electric 
(TE)  (magnetic  field  in  the  z  direction).  These  two  types  of 
modes  can  be  described  by  two  scalar  wave  equations  [14]. 
Assuming  periodicity  of  the  solution  in  the  xy  plane,  the 
Floquet-Bloch  theory  shows  that  the  scalar  fields  satisfy  a 
Hermitian  eigenvalue  equation  in  a  unit  cell  of  the  form 

Au—XAAu  in  £2.  (1) 

For  the  TE  mode,  we  have  u  =  Hz(r ),  X  =  co2 /c2,  and 

A(e,k)  =  -(V  +  ik)  •  [e_1(r)(V  +  ik)],  M  =  I, 

where  c  is  the  speed  of  light,  I  denotes  the  identity  operator, 
and  the  wave  vector  k  lies  in  the  irreducible  Brillouin  zone  B. 
For  the  TM  mode,  we  have  u  =  Ez(r ),  X  =  co2 /c2,  and 

Aik)  =  -(V  +  ik)  •  (V  +  ik),  M(e)  =  e(r)I. 

The  unit  cell  Q  and  the  irreducible  Brillouin  zone  B  depend 
on  the  lattice  type  [5].  We  denote  by  ( um,Xm )  the  mth 
eigenfunction-eigenvalue  pair  of  Eq.  (1)  and  assume  that  these 
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eigenpairs  are  numbered  in  ascending  order:  0  <  X1  <  X2  < 
...^A°°. 

Let  J  —  {nij  |  j  =  1,2 ,...,/}  denote  a  set  of  J  bands  for 
which  we  seek  to  achieve  absolute  gaps.  For  instance,  J  = 
{1,3,5}  indicates  the  set  of  the  first,  third,  and  fifth  band  gaps 
corresponding  to  mi  =  1, m2  =  3,  and  m3  =  5.  We  define  the 
eigenvalue  gap-midgap  ratio  between  bands  m7  and  m7  +  1 
as 

inf  Xmj+l(€(r),k)  -  sup Xm^(e(r),k) 

keB  IcgB 

RJ(e(r))  =  2 - - - — - . 

inf  Xmj+l(e(r),k)  +  sup  Xmj(e{r),  k) 

keB  keB 

We  therefore  consider  the  following  optimization  problem: 

sup  inf  a;  Rj(€(r))  .  (2) 

Here  Q  =  |6(r)  :  e(r)  e  [6^,6#],  Vr  e  £2}  is  the  admissible 
domain,  where  eL  and  are  dielectric  constants  of  a  low- 
index  material  and  a  high-index  material,  respectively.  The 
cij  are  prescribed  weights  for  each  band  gap  (though  in  our 
computation  we  typically  set  otj  =  1,  j  =  1, . . . ,/).  Thus  the 
objective  in  Eq.  (2)  is  to  find  an  optimal  material  distribution 
that  maximizes  the  (weighted)  smallest  gap  among  J  chosen 
bands. 

In  practice,  we  discretize  the  optimization  problem  (2)  as 
follows.  First,  we  consider  only  n k  wave  vectors  in  the  set 

Vh  =  {kk  e  dB,  1  <  k  ^  nk },  (3) 

where  3  B  represents  the  boundary  of  the  irreducible  Brillouin 
zone.  Second,  the  unit  cell  Q  is  discretized  into  N  x  N 
elements  on  which  the  dielectric  function  takes  a  piecewise 
constant  value  between  eL  and  6#  on  each  element;  moreover, 
if  the  symmetry  of  the  prescribed  lattice  is  taken  into 
consideration,  the  dielectric  function  only  needs  to  be  defined 
in  1/8  of  the  unit  cell  in  a  square  lattice,  or  in  1/12  of  the  unit 
cell  in  a  hexagonal  lattice,  namely, 

Q^{6:6G[6L,6^},  (4) 

where  n€  <  N2  =  J\f.  Third,  we  use  a  Galerkin  finite  element 
method  with  piecewise  linear  polynomials  to  approximate  the 
system  (1)  as 

Ah(e,k)umh]  =  Xmhj Mh{e)umhj ,  6  e  Qh,k  e  Vh,  (5) 

where  Ah{e,k)  e  C^xAf  is  a  Hermitian  matrix  and  Mh(e )  e 

RWxV 

is  a  symmetric  positive  definite  matrix.  Since  e(r)  is 
piecewise  constant  on  Q,  the  6 -dependent  matrices  can  be 
expressed  as 

A™(t,k)  =  jr  €-'A^(k),  M™(()  =  jr  etM™  (6) 

i— 1  i= 1 

while  A™(k)  and  M^E  are  independent  of  6.  Here  and  below, 
superscripts  ™  and  TE  are  used  to  indicate  TM  polarization 
and  TE  polarization,  respectively. 
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DESIGN  OF  PHOTONIC  CRYSTALS  WITH  MULTIPLE  AND  . . . 


For  the  TE  case,  we  define  the  discrete  eigenvalue  gap- 
midgap  ratio  between  eigenvalues  m  j  and  Mj  +  1  as 


dT  EJ 

Kh 


(0  =  2 


min  X 

keVh 


TE,m7  +  l 
h 


(e,  k)  —  max  Xh  ’  J 

keVh  n 


(c,k) 


min  A. 

keVh 


TE,m;  +  l 
h 


(e,k)  +  max  xTE,mj 
keVh 


(c,k) 


For  the  TM  case  the  discrete  eigenvalue  gap-midgap  ratio 
R™J(e )  is  defined  in  a  similar  manner.  For  the  complete  band 
gaps  we  define  the  discrete  eigenvalue  gap  ratio  between  bands 
mj  and  mj  +  1  as 


RJh(0  =  min  [Rh  /?,  '(>)]. 

To  design  the  photonic  crystal  structure  that  supports  multiple 
combined  band  gaps  we  propose  to  solve  the  following 
optimization  problem: 


max  min  otj  RJh{ e),  s.t. 


A^{e,k)u 


TE, mj 


=  X 


TE  ,m , 


MT(e)uh 


TE, mj 


A™(e,k)u™’mj  =  X™’mj  M™(e)u 
for  j  =  l,...  ,J,k  e  Vh- 


TM,? 

h 


(7) 


To  model  only  TE  or  TM  polarization,  one  can  simply  omit  the 
nonrelevant  equations  in  the  above  formulation.  Note  further 
that  the  formulation  (7)  can  be  generalized  to  treat  more 
general  cases  in  which  the  number  and  location  of  TE  band 
gaps  are  allowed  to  differ  from  those  of  TM  band  gaps. 

In  the  previous  work  [14],  we  showed  how  to  use  semidef- 
inite  inclusions  combined  with  subspace  methods  to  locally 
approximate  the  single  ABG  problem  [only  one  band  gap  and 
either  TE  or  TM  polarization,  which  is  simpler  than  Eq.  (7)] 
as  a  convex  semidefinite  program;  see  Sec.  3.3.1  of  [14].  We 
now  indicate  how  extensions  of  these  ideas  can  be  used  to 
develop  a  tractable  (i.e.,  conic  and  convex)  local  approximation 
of  Eq.  (7).  It  will  be  convenient  to  consider  the  change  of 
variables  defined  by  yj  =  l/e7- ,  j  =  1 ,  ... ,n€.  Let  £  be  a  given 
parameter  vector  satisfying  £  e  Qh  and  define  yj  =  1/ey, 
j  =  l,  ...  ,n€.  Then  the  following  optimization  problem  is 
an  extension  of  approximations  (16)  and  (17)  developed 
in  [14]: 

max  F 

i  ^*(Y,k)[AT(Y,k)  -  bjM™]<t>™(y,k)  <  0, 
*™*(y,k)[AlE(y,k)  -  ajtfF]9™(9,k)  h  0, 
4>™*(€,k)[djA™(k)  -  ,k)  <  0, 

'i >™*(Ck)[cjA™(k)  -  M™(()]*TmM((,k)  >  0, 


CiYi  =  1, 

aj  ^  0,  bj  ^  0,  Cj  ^  0,  dj  ^  0, 

e  e  Qh,y  e  Sh,k  e  Vh,  and  j  e  {1, . . . ,/}, 
i  =  1,  ...  ,n€, 

where  y  =  \e,y,a,b,c, d,F]  and  Sh  =  {y  :  y  G 


As  developed  in  [14],  the  matrices  <£^(j> ,  k)  and  ty^(y,k), 
[respectively,  $>™(e,k)  and  ^™(€,k)]  ideally  are  comprised 
columnwise  of  the  lower  mj  eigenvectors  and  the  upper 
M  —  mj  eigenvectors  of  Eq.  (5)  for  the  TE  case  (respectively, 
the  TM  case).  As  developed  in  [14],  we  instead  work  with 
a  small  “important”  subset  of  these  eigenvectors  to  keep  the 
computation  efficient;  see  [14]  for  details.  The  semidefinite 
inclusions  [the  first  four  sets  of  constraints  in  Eq.  (8)]  are  linear 
in  the  design  variables  y.  The  bilinear  constraints  [the  fifth 
through  seventh  set  of  constraints  in  Eq.  (8)]  can  be  linearized 
around  the  previous  solution  to  obtain  a  linear  semidefinite 
program.  For  instance,  the  constraint  e*  yi  =  1  is  linearized  as 
fi  +  €t(yi  —  ft)  +  fife  —  €i)  =  1.  The  resulting  linearized 
semidefinite  program  can  be  efficiently  solved  using  modern 
interior  point  methods  such  as  the  SDPT3  software  [16]. 


III.  MESH  ADAPTIVITY  AND  COMPUTATIONAL 
PROCEDURE 

We  develop  an  adaptive  mesh  refinement  algorithm  in 
order  to  improve  the  efficiency  and  accuracy  of  discretization. 
Rather  than  using  a  uniform  computational  mesh,  we  employ 
a  discretization  adapted  to  the  particular  material  distribution 
that  is  being  computed.  Given  any  initial  coarse  representation 
of  the  optimal  configuration,  we  choose  and  subdivide  ele¬ 
ments  that  meet  the  refinement  criteria.  In  particular,  elements 
are  refined  when  they  are  on  the  material  interface,  i.e., 
when  there  is  a  jump  in  the  value  of  the  design  variables 
between  the  element  under  consideration  and  at  least  one  of 
its  neighbors.  When  an  element  is  refined,  hanging  nodes 
are  generated  if  the  neighboring  elements  are  of  different 
sizes  after  refinement.  In  order  to  simplify  the  computational 
procedure,  we  impose  a  2: 1  rule  to  restrict  the  refinement  level 
difference  between  neighboring  elements:  if  the  refinement 
level  difference  exceeds  2,  the  larger  element  is  further  refined 
such  that  at  most  one  hanging  node  exists  on  any  element  edge. 
In  order  to  ensure  a  continuous  finite  element  interpolation, 
the  degree  of  freedom  corresponding  to  the  hanging  node  is 
constrained  to  the  interpolated  value  from  the  two  comer  nodes 
and  can  therefore  be  locally  eliminated  from  the  global  system. 

Our  computational  procedure  incorporating  mesh  adaptiv¬ 
ity  in  conjunction  with  the  local  optimization  problem  (8)  is 
as  follows: 

(1)  Start  with  a  coarse  mesh  8x8  and  a  random  initial 
distribution  £,  and  an  error  tolerance  etoi; 

(2)  if  necessary,  use  interpolation  to  adapt  the  current 
distribution  e  to  correspond  to  the  current  refined  mesh; 

(3)  for  each  wave  vector  k  e  Ph,  compute  the  subspaces 
<&™(i,k),  (ifk),  0™(y,k),  and  ty™(y,k),  and  form  the 
linearized  version  of  the  semidefinite  program  (8); 

(4)  solve  the  linearized  version  of  Eq.  (8)  for  an  optimal 
solution  e*; 

(5)  if  ||  6*  —  i  ||  ^  6toi,  and  current  refinement  level  < 
maximum  refinement  level ,  then  refine  elements  adaptively, 
extrapolate  £  of  the  new  mesh,  and  go  to  1; 

(6)  else  if  \\e*  —  £\\  >  €to\ ,  and  current  refinement  level  ^ 

maximum  refinement  level ,  then  update  £  €*,  and  go  to  2; 

(7)  else  if  \\e*  —  £\\  ^  etoi,  and  current  refinement  level  ^ 
maximum  refinement  level ,  then  stop  and  return  optimal  €*. 
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IV.  RESULTS  AND  DISCUSSION 

We  present  computational  results  for  the  design  of  photonic 
crystals  made  up  of  doped  GaAs  (er  =  12.8)  and  vacuum 
air  (€a  =  1)  [10,15].  We  computed  a  wide  variety  of  opti¬ 
mized  pairs  and  triplets  of  multiple  band  gaps  in  TE,  TM, 
and  complete  TEM  modes,  for  both  square  and  hexagonal 
lattices.  In  the  results  presented  herein,  the  eigenvalues  are 
plotted  in  the  dimensionless  unit  (coa/lnc)2,  where  a  is  the 
lattice  constant  for  both  square  and  hexagonal  lattices.  For 
comparison  between  our  results  and  other  previous  results,  we 
also  describe  our  results  in  terms  of  th q  frequency  gap-midgap 
ratio  between  bands  mj  and  mj  +  1,  defined  as 

A  mj  min  co?J+1(€,k)  —  max  ox  7  (c, A:) 

nj  =  AV  k£Pn  h_  kJn  h_ 

—  —mj  Z  .  m.  +  l  mj  • 

°>h  min  coh  (€,k)  +  max  ox  (e,k) 

keVh  n  keVh 

The  frequency  gap-midgap  ratio  is  used  as  the  objective 
function  in  some  published  research  (see,  e.g.,  [10,15]),  while 
the  eigenvalue  gap-midgap  ratio  is  used  in  other  published 
research  (e.g.,  [9,13,14]).  Despite  the  obvious  difference 
between  the  eigenvalue  and  frequency  (the  former  being  the 
square  of  the  latter  divided  by  the  speed  of  light),  the  optimal 
crystal  structures  have  been  observed  to  be  astonishingly 
consistent  when  either  is  used  in  the  gap  objective  function  in 
single  band  gap  optimization  problems  [10,13,14].  While  it  is 
intuitive  that  the  frequency  relative  gap-midgap  ratio  should  be 
monotone  in  the  eigenvalue  gap-midgap  ratio,  one  can  create 
pathological  counterexamples.  Herein  we  choose  to  optimize 
the  eigenvalue  gap-midgap  ratio  because  the  first  four  sets 
of  constraints  of  our  optimization  model  (8)  are  linear  in  the 
eigenvalues  and  so  require  no  extra  linearization  themselves, 
and  the  fifth  and  sixth  constraints  are  only  modestly  nonlinear 
in  the  eigenvalues.  Of  course,  should  one  wish  to  optimize  the 
frequency  gap-midgap  ratio,  the  resulting  nonlinear  constraints 
[the  fifth  and  sixth  set  of  constraints  in  Eq.  (8)]  could  be 
linearized  as  discussed  earlier  when  constructing  the  linear 
semidefinite  program. 

We  initialize  our  procedure  with  a  very  coarse  grid  8x8 
and  adaptively  refine  it  up  to  128  x  128  in  resolution.  We 
initialize  with  a  random  distribution  e  on  this  coarse  grid, 
and  then  obtain  an  optimal  solution  on  the  coarse  grid  8x8; 
this  optimal  solution  is  used  as  the  basis  to  refine  the  grid 
and  is  then  adapted  to  the  refined  grid  by  interpolation  as 
discussed  in  Sec.  III.  (A  detailed  discussion  on  the  choice  of 
some  other  simulation  parameters  can  be  found  in  our  previous 
work  [14]).  We  found  that  the  grid  refinement  procedure 
significantly  reduces  both  the  number  of  degrees  of  freedom 
in  the  finite  element  procedure  and  the  number  of  decision 
variables  in  the  optimization  procedure,  thereby  leading  to  a 
significant  reduction  in  computation  time,  with  a  speedup  in 
overall  computation  in  the  range  40%  -500%,  but  typically 
around  250%.  All  computations  were  performed  on  a  Linux 
PC  with  Dual  Core  AMD  Opteron  270,  2.0  GHz,  and  each  run 
of  our  procedure  was  obtained  in  1-10  min.  The  relatively  low 
computation  times  enabled  us  to  study  the  inherent  tradeoffs 
between  optimizing  two  different  band  gaps  (namely,  the 
tradeoff  frontier),  which  provides  very  useful  information 
for  choosing  the  most  appropriate  design.  Below  we  present 
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representative  results  to  illustrate  various  aspects  of  solutions 
to  the  multiple  and  combined  band  gap  optimization  problem. 

A.  TM  band  gaps 

We  illustrate  a  typical  sequence  of  optimal  solutions  using 
our  mesh  adaptivity  procedure  on  the  problem  of  optimizing 
the  second  and  fifth  band  gaps  [with  identical  weights,  i.e., 
(q'1,0'2)  =  (1,1)]  in  the  hexagonal  lattice.  Our  results  are 
shown  in  Fig.  1.  It  is  important  to  note  that  grid  resolution 
affects  the  obtained  results.  In  particular,  relatively  coarse 
resolutions  ( h min  =  <2/8, a /l 6)  produce  optimal  solutions  with 
some  mixed  features  (i.e.,  €l  <  G  <  £h  for  some  cell  elements 
/),  while  finer  resolutions  (hm  in  =  a/32,  a/64,  a/128)  yield 
optimal  structures  which  involve  only  pure  concentrations 
(6/  =  eL  or  eH  for  all  cell  elements  /).  Moreover,  the  shape 
of  the  circular  inclusions  is  much  more  visible  on  the  finer 


FIG.  1.  Mesh  adaptivity  results  show  the  computational  grids 
overlaid  on  the  crystal  structures  (left),  optimal  crystal  structures 
(middle),  and  frequency  bands  (right)  for  the  second  and  fifth  TM 
band  gaps  in  the  hexagonal  lattice.  Grid  resolution  varies  from 
hrmn  =  a/ 8  (top  panel),  /imin  =  a/ 16  (second  panel),  /imin  =  a/32 
(third  panel),  h ^  =  a/64  (four  panel),  to  hm in  =  a/ 128  (bottom 
panel).  The  letters  T,  M,  and  K  in  this  figure  and  some  figures 
below  represent  the  vertices  of  the  irreducible  Brillouin  zone  of  the 
hexagonal  lattice. 
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FIG.  2.  (Color  online)  Optimization  results  show  the  final  compu¬ 
tational  grid  (left),  optimal  crystal  structure  (middle),  and  frequency 
bands  (right)  for  the  first,  second,  and  fourth  TM  band  gaps  in  the 
square  lattice.  The  letters  T,  X,  and  M  in  this  figure  and  some  figures 
below  represent  the  vertices  of  the  irreducible  Brillouin  zone  of  the 
square  lattice. 


grids.  Our  optimized  structure  differs  from  existing  single-gap 
photonic  crystals  [10,13,14]  in  the  sense  that  it  has  both 
small  inclusions  (of  radius  rs/a  =  0.06)  and  large  inclusions 
(of  radius  ri/a  =  0.12),  while  single-gap  photonic  crystals 
typically  consist  of  inclusions  of  the  same  radius.  It  is 
interesting  to  observe  that  the  unit  cell  of  the  optimized 
structure  consists  of  two  large  circular  disks  and  three  small 
disks  as  we  optimize  the  second  and  fifth  band  gaps. 

A  typical  triplet  of  optimized  band  gaps  is  shown  in  Fig.  2, 
which  presents  the  optimized  structures  for  optimizing  the  first, 
second,  and  fourth  band  gaps  (with  identical  weights  ot\  = 
a2  =  a?3  =  1)  in  the  square  lattice.  In  this  case,  the  optimized 
structure  also  consists  of  circular  disks  of  different  sizes 
{ri/a  =  0.17 ,rs/a  =  0.08).  Similar  to  the  previous  results, 
the  relative  eigenvalue  gaps  are  the  same  for  all  three  bands 
since  we  chose  the  same  weights  for  all  bands.  One  notable 
feature  of  this  photonic  structure  is  that  the  midgap  frequency 
of  the  third  band  is  approximately  twice  that  of  the  first  band. 

In  optimizing  a  weighted  pair  of  band  gaps,  there  is  an 
intuitive  tradeoff  between  the  size  of  one  band  gap  versus  the 
other  band  gap  that  can  be  computed  by  varying  the  weights 
associated  with  each  band  to  yield  a  tradeoff  frontier  between 
the  two  band  gaps.  Such  a  tradeoff  frontier  is  illustrated  in 
Fig.  3  for  the  problem  of  optimizing  the  weighted  first  and  third 


FIG.  3.  The  tradeoff  frontier  for  the  first  (horizontal  axis)  and 
third  (vertical  axis)  TM  band  gaps  in  the  hexagonal  lattice.  The 
frontiers  are  traced  by  varying  the  weights  corresponding  to  each  band 
gap.  Multiple  frontiers  can  be  attributed  to  multiple  local  optima. 


band  gaps  in  the  hexagonal  lattice.  The  points  in  the  figure 
were  produced  by  varying  the  weights  (a\,a2)  associated 
with  two  band  gaps  (first  and  third  bands,  respectively)  for 
a  variety  of  values  of  (ai,a2)  e  [0,  l]2.  More  specifically,  we 
start  by  computing  a  solution  for  (ai,a2)  =  (0.5, 0.5)  on  a 
uniform  mesh  64  x  64.  We  then  modify  the  weights  (a\,a2)  e 
[0,  l]2  slightly  and  use  the  previously  obtained  solution  as 
the  initiating  distribution  for  computing  the  solution  for  the 
problem  with  modified  weights.  This  step  is  repeated  in  order 
to  track  various  local  optimal  branches.  Figure  3  illustrates  the 
tradeoff  frontiers  A -B  and  C-D,  where  A,  B ,  C,  and  D  cor¬ 
respond  to  the  optimized  structures  for  (oq  ,a2)  =  (0.99,0.01), 
(0.01,0.99),  (0.99,0.01),  and  (0.01,0.99),  respectively.  The 
figure  shows  two  different  frontiers  A -B  and  C-D.  The 
multiple  branches  of  the  frontiers  represent  multiple  significant 
local  optima  corresponding  to  identical  weights.  Structure  A 
favors  the  third  band  gap,  whereas  structure  B  favors  the  first 
band  gap.  Structure  C  exhibits  a  good  compromise  between 
the  two  band  gaps  since  both  are  relatively  large  in  this  case. 
Structure  C  is  particularly  interesting  in  that  it  resembles  D  in 
overall  design,  but  has  larger  disks  with  an  air  hole.  Frontier 
A -B  consists  of  structures  having  disks  of  two  different  radii 
but  otherwise  similar  topology,  while  frontier  C-D  consists  of 
structures  having  disks  of  different  radii  and  different  topology. 


B.  TE  band  gaps 

We  now  turn  to  TE  band  gaps.  Figures  4  and  5  show  typical 
results  for  optimizing  pairs  (first  and  second  band  gaps)  and 
triplets  (second,  fourth,  and  sixth  band  gaps)  in  the  square 
lattice,  respectively,  where  all  bands  have  equal  weights.  We 
observe  that  the  optimized  structure  in  Fig.  4  is  connected  and 
relatively  simple,  whereas  the  optimized  structure  in  Fig.  5 
is  nonconnected  and  far  more  complicated.  Although  TE 
polarization  typically  favors  connected  lattices,  these  results 
show  that  it  is  possible  to  obtain  nonconnected  TE  structures 
with  multiple  band  gaps. 

Figure  6  shows  solutions  for  optimizing  the  first  and  fourth 
TE  band  gaps  for  both  the  square  and  hexagonal  lattices 
(with  equal  weights  for  the  two  bands).  We  observe  that  the 
optimized  structure  in  the  square  lattice  is  not  connected,  while 
the  optimized  structure  in  the  hexagonal  lattice  is  connected. 
Moreover,  we  see  for  both  band  structures  that  the  midgap 
frequency  of  the  fourth  band  is  roughly  twice  that  of  the  first 
band.  Therefore  similar  to  the  structure  shown  in  Fig.  2,  these 
designs  can  prohibit  electromagnetic  waves  at  both  frequencies 
co  and  2 co. 

We  also  studied  the  tradeoff  frontier  for  the  first  and 
third  TE  band  gap  in  the  hexagonal  lattice;  see  Fig.  7.  This 


FIG.  4.  Optimization  results  show  the  final  computational  grid 
(left),  optimal  crystal  structure  (middle),  and  frequency  bands  (right) 
for  the  first  and  second  TE  band  gaps  in  the  square  lattice. 
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FIG.  5.  Optimization  results  show  the  final  computational  grid 
(left),  optimal  crystal  structure  (middle),  and  frequency  bands  (right) 
for  the  second,  fourth,  and  sixth  TE  band  gaps  in  the  square  lattice. 


frontier  has  proven  to  be  more  complicated  than  its  TM 
mode  counterpart  shown  previously  in  Fig.  3.  In  fact,  no 
distinctive  frontier  can  be  observed,  which  is  undoubtedly 
due  to  numerous  local  optima  in  this  case.  We  employed  a 
similar  computational  strategy  to  that  used  to  produce  Fig.  3  as 
described  in  the  previous  subsection.  In  Fig.  7  we  display  four 
(locally)  optimized  structures  along  the  envelope  A-B-C-D 
of  the  apparent  frontier,  where  A,  B ,  C,  and  D  correspond  to 
solutions  for  (ai,a2)  =  (0.99,0.01),  (0.5, 0.5),  (0.5, 0.5),  and 
(0.01,0.99),  respectively.  Structure  A  favors  the  third  band, 
whereas  structure  D  favors  the  first  band.  Structures  B  and 
C  offer  a  compromise  between  the  two  bands.  Note  further 
that  the  two  structures  B  and  C  are  very  different  despite  the 
fact  that  they  have  similar  objective  values.  (We  also  observed 
many  local  optimal  solutions  in  our  previous  work  on  the 
optimal  design  of  photonic  crystals  with  single  band  gap  [14].) 


C.  Complete  band  gaps 

We  are  able  to  compute  complete  band  gaps  for  both  the 
hexagonal  and  square  lattices;  the  band  gaps  in  the  hexagonal 
lattice  lie  in  the  first  and  third  bands  for  the  TE  polarization, 
yet  in  the  third  and  sixth  bands  for  the  TM  polarization. 
Figure  8  illustrates  the  geometry  and  band  structures  for  the 


FIG.  6.  Optimization  results  show  optimal  crystal  structure  (left), 
and  frequency  bands  (right)  for  the  first  and  fourth  TE  band  gaps  in 
the  square  lattice  (top)  and  the  hexagonal  lattice  (bottom). 
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FIG.  7.  The  tradeoff  frontier  for  the  first  (horizontal  axis)  and 
third  (vertical  axis)  TE  band  gaps  in  the  hexagonal  lattice.  The 
frontiers  are  traced  by  varying  the  weights  corresponding  to  each  band 
gap.  Multiple  frontiers  can  be  attributed  to  multiple  local  optima. 

hexagonal  lattice.  The  corresponding  frequency  gap-midgap 
ratios  of  5.76%  and  6.94%  are  quite  sizable — and  are  the  first 
multiple  and  complete  band  gaps  ever  reported  for  photonic 
crystals  in  the  hexagonal  lattice.  This  photonic  crystal  structure 
has  both  connected  and  nonconnected  features.  Moreover,  it 
has  a  more  complicated  geometry  than  the  previous  structures 
shown  herein. 

For  the  square  lattice,  complete  band  gaps  lie  in  the  third 
band  for  the  TE  polarization,  yet  in  the  sixth  and  ninth  bands 
for  the  TM  polarization.  Figure  9  illustrates  our  results.  The 
corresponding  frequency  gap -midgap  ratios  of  7.59%  and 
13.5%  are  also  large.  We  emphasize  again  that  no  other 
multiple  and  complete  gaps  have  previously  been  found  for 
photonic  crystals  in  the  square  lattice.  In  general,  the  complete 
band  gaps  are  smaller  than  either  of  the  corresponding  TE  and 
TM  band  gaps  because  it  is  rather  difficult  to  simultaneously 
achieve  both  the  TE  and  TM  band  gaps.  Of  course,  one  can  also 
widen  the  gap  size  for  one  band  at  the  expense  of  narrowing 
the  gap  size  of  the  other  band  by  choosing  appropriate  weights 


FIG.  8.  Optimization  results  for  the  multiple  complete  band  gaps 
in  the  hexagonal  lattice.  The  optimal  crystal  structure  is  shown  in  the 
left-hand  picture;  while  the  frequency  bands  are  shown  in  the  right- 
hand  picture,  with  solid  lines  representing  TM  bands,  and  broken 
lines  representing  TE  bands.  The  first  complete  band  gap  is  formed 
by  the  overlapping  of  the  first  TE  and  third  TM  band  gaps,  while  the 
second  complete  band  gap  is  formed  by  the  third  TE  and  sixth  TM 
band  gaps. 
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FIG.  9.  Optimization  results  for  the  multiple  complete  band  gaps 
in  the  square  lattice.  The  optimal  crystal  structure  is  shown  in  the 
left-hand  picture;  while  the  frequency  bands  are  shown  in  the  right- 
hand  picture,  with  solid  lines  representing  TM  bands,  and  broken 
lines  representing  TE  bands.  The  first  complete  band  gap  is  formed 
by  the  overlapping  of  the  third  TE  and  sixth  TM  band  gaps,  while  the 
second  complete  band  gap  is  formed  by  the  third  TE  and  ninth  TM 
band  gaps. 

in  the  optimization  problem  (8).  Finally,  it  is  interesting  to  note 
in  this  case  that  although  the  photonic  crystal  structure  has  a 
complicated  geometry,  it  is  nevertheless  connected. 

V.  CONCLUSIONS 

In  conclusion,  we  have  demonstrated  the  usefulness  of  a 
computational  scheme  based  on  conic  (semidefinite)  convex 
optimization  to  design  photonic  crystals  with  multiple  and 
complete  band  gaps  in  both  square  and  hexagonal  lattices. 
Our  photonic  crystal  patterns  are  different  from  existing 
photonic  crystal  designs  [10,13,14]  with  a  single  band  gap. 


In  particular,  we  observe  from  our  results  that  unlike  optimal 
photonic  crystals  discovered  in  [13],  photonic  crystals  with 
multiple  band  gaps  do  not  follow  some  simple  geometric 
properties.  Therefore  numerical  optimization  is  crucial  to 
the  design  of  photonic  crystals  that  support  several  band 
gaps.  These  photonic  crystals  may  prove  useful  for  the 
suppression  of  resonance  at  harmonic  frequencies,  as  they 
prohibit  the  propagation  of  electromagnetic  waves  at  several 
different  frequencies.  In  addition,  we  have  computed  photonic 
structures  with  large  complete  band  gaps.  Our  results  hopefully 
open  up  a  new  arena  for  the  design  of  photonic  crystals. 

We  note  that  many  of  the  optimized  crystal  designs 
shown  herein  involve  intricate  patterns  of  materials  at  the 
nano-level,  and  may  be  too  expensive  or  even  impossible 
to  fabricate.  Simply  incorporating  “fabrication  constraints” 
such  as  connectedness  of  materials  or  bounds  on  the  curva¬ 
ture  of  boundaries  easily  yields  combinatorially  intractable 
optimization  models.  Instead,  we  propose  to  modify  the  basic 
optimization  problem  so  that  a  resulting  solution  is  “robust” 
for  fabrication,  somewhat  in  the  spirit  of  robust  convex 
optimization  [17].  This  is  the  subject  of  ongoing  research. 
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